Identifying reactive trajectories using a moving transition state 
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A time-dependent no-recrossing dividing surface is shown to lead to a new criterion for identifying 
reactive trajectories well before they are evolved to infinite time. Numerical dynamics simulations of 
a dissipative anharmonic two-dimensional system confirm the efficiency of this approach. The results 
are compared to the standard fixed transition state dividing surface that is well-known to suffer from 
recrossings and therefore requires trajectories to be evolved over a long time interval before they 
can reliably be classified as reactive or non-reactive. The moving dividing surface can be used to 
identify reactive trajectories in harmonic or moderately anharmonic systems with considerably lower 
numerical effort or even without any simulation at all. 



I. INTRODUCTION 

It is perhaps surprising how many problems in chem- 
istry, physics, and biology can be reduced to the simple 
model of diffusion over a barrier .Ij . Although chemical 
reactions in all phases of matter provide the prime exam- 
ple [1, 131 J a plethora of systems that evolve from suitably 
defined "reactant" to "product" states are amenable to 
a description in this framework |1, IE IE Q, Hi Q- Tran- 
sition State TheorWTST) H Efl |nj or one of its de- 
scendants [l3, Us nil is often used to approximate the 
rate of these reaction processes. This theory is based on 
the assumption that the reaction rate is determined in 
a small volume of the phase space near the barrier. It 
is then possible to define a dividing surface separating 
reactants from products and obtain the rate from the 
flux through this surface. The optimal location of the 
dividing surface is that which minimizes the number of 
recrossings — the fundamental idea of variational transi- 
tion state theory EHIl E El El . When the system of 
interest can be viewed as isolated from its environment, 
as in low-density gas phase chemical reactions, TST may 
indeed provide an excellent approximation to the rate. 
However, most processes of interest do not take place in 
isolation, but rather in a complex environment where in- 
teractions between the system and its surroundings occur 
on time scales comparable to that of the reaction: In a 
reaction occurring in the condensed phase, for instance, 
the dynamics of the solute is strongly coupled to that of 
the solvent. In this case, the fundamental assumption of 
TST, that the dividingsurface is crossed once and only 
once, no longer holds ;23,|2l|,|22,|23|- On the time scale of 
the reaction event, fluctuations of the environment will 
almost inevitably cause recrossings of the dividing sur- 
face that lead to an overestimation of the rate. 

To overcome the recrossing problem, many TS di- 
viding surfaces have been suggested in the literature 
[3, 19, 24, 25j which provide systematic (and simple) 
approximations to the optimal TS dividing surface. In 
some cases, the dividing surface has been identified in 



the infinite-dimensional phase space consisting of the sys - 
tem and an explicit set of bath oscillators |2E |2a, |23| ■ 
This approach leads to an excellent approximation to 
the rate [Tj, |23 • Interestingly, the same result was sub- 
sequently obtained without recourse to the explicit heat 
bath model, using instead a collective reaction coordinate 
containing the influence of the bath directly [23 ■ 

In a recent series of papers |23, Is^ , we reformulated 
the recrossing problem using a dividing surface that is 
itself moving stochastically so as to avoid recrossings. 
The motion of that surface follows the unique trajec- 
tory — named the TS trajectory — that never leaves the 
barrier region. Any reactive trajectory crosses the mov- 
ing surface once and only once, whereas a nonreactive 
trajectory does not cross at all. This construction ex- 
tends the approach of [23 in that it provides not only 
a reaction coordinate, but also the complete geometric 
structure by specifying all of the unstable and stable 
degrees of freedom glob ally [31|. The previous purely 
analytic studies |29l l30l | are complemented here with a 
numerical investigation of the reaction dynamics for a 
two-dimensional stochastic nonlinear model. It will be 
shown that the moving dividing surface offers consider- 
able computational advantages over the traditional fixed 
surface: Its use can significantly reduce the simulation 
time required to distinguish between reactive and non- 
reactive trajectories. Indeed, for a harmonic barrier it 
identifies reactive trajectories a priori, so that the need 
to simulate their dynamics does not arise at all. In an 
anharmonic system, the identification of reactive trajec- 
tories by the moving surface is no longer exact. Neverthe- 
less, for moderately strong anharmonicities it provides a 
useful approximation, and its advantages over the fixed 
surface are retained. In addition, the moving TS surface 
introduces novel observables that characterize the reac- 
tion process on a microscopic level. Most prominently, 
it allows one to define a unique reaction time for each 
individual trajectory. 

The outline of this paper is as follows: In Sec. m the 
two-dimensional stochastic nonlinear model that is the 



focus of the computational studies in this work is de- 
fined and the construction of the moving TS dividing 
surface and its associated geometric structures is briefly 
reviewed. In Sec. lIIII the ensemble of trajectories is spec- 
ified by a thermal distribution of particles localized at the 
conventional TST dividing surface. This barrier ensem- 
ble is reminiscent of the weighting distribution in stan- 
dard rate expressions and is appropriate even in nonlinear 
cases. Its simple structure also readily leads to the an- 
alytic determination of several observables of the model 
system in the harmonic limit (Sec lIVJI . They are in pre- 
cise agreement with the numerical results presented in 
Sec. The latter section also demonstrates that ob- 
servables converge faster when evaluated using the mov- 
ing dividing surface rather than conventional numerical 
methods, both in the harmonic limit and in systems with 
anharmonic barriers. This observation is particularly 
useful in the anharmonic case when the chosen system 
is not amenable to analytic approaches. 



II. PRELIMINARIES 

Although the general theory is applicable to systems 
with an arbitrary number n of degrees of freedom, the fol- 
lowing discussion will be restricted to n = 2 coordinates 
under the influence of a stochastic bath. This choice can 
be made without loss of generality because it exhibits 
all the salient features of the higher-dimensional cases: 
It can encompass an unbound (reactive) direction and 
a bound bath mode whose interaction with the reactive 
mode is strong enough to require its explicit description. 
The coupling of the modes is described by a Taylor ex- 
pansion about a transition point (or col) on the potential 
energy surface. Such a model with a minimum number 
of nonlinear terms is described in this section. It will be 
used in the following to study the effect of increasing an- 
harmonicity on the identification of reactive trajectories. 

To set the stage for the following investigations, the 
construction of the moving TST dividing surface and the 
associated invariant manifolds is summarized in the re- 
mainder of this section. The reader interested in a full 
exposition is referred to Refs. |23 andlSOl where the for- 
malism was first introduced. 



A. The Tw^o-Dimensional Dissipative Model 

A prototypical reactive system within a solvent may 
be described by the Langevin equation [SjI 



force. The latter is related to the friction matrix T by 
the fluctuation-dissipation theorem [33 



ga(i) = -V,-C/(g;(i)) - Tq^it) + ^„(i) 



(1) 



The vector q here denotes a set of n = 2 mass- weighted 
coordinates, U{q} the potential of mean force governing 
the reaction, T a symmetric positive-deflnite friction ma- 
trix and Cait) a fluctuating force assumed to be Gaus- 
sian with zero mean. The subscript a represents ran- 
domness by labeling different instances of the fluctuating 



Ut)^it')) =2kBTTS{t-t') 



(2) 



where the angular brackets denote the average over the 
instances a of the noise. Although not strictly necessary, 
the friction is often taken to be isotropic, i.e., 



r = 7J 



(3) 



with a scalar friction constant 7. 

The reactant and product regions in configuration 
space are separated by a potential barrier whose posi- 
tion is marked by a saddle point % = oi the potential 
U{q). In its vicinity, the potential is approximately har- 
monic and can always be written in a diagonal normal 
form. In general, anharmonic terms will be present in 
the potential. In the neighborhood of the saddle point, 
where the reaction rate is determined, they are only mod- 
erately strong, but usually not negligible. In this work, 
we include a typical (even) higher-order nonlinear term 
and focus on the potential 



U{x, y) = --ujIx^ + -ujly' + kx^y" 



(4) 



where the position vector is written as q — {x,y), and 
the constant k quantifies the nonlinear coupling of the 
different degrees of freedom. Note that the nonlinear- 
ity in the potential (0)) is symmetric in the coordinate 
system and neglects other fourth order terms that are 
typically retained in the analysis of anharmonic barriers. 
(See, e.g., Ref.l33. in which such coupled anharmonic po- 
tentials have been used to study the H -I- II2 ^ H2 -I- H 
reaction and bound vibrational systems.) However, as 
discussed in the Appendix, it is amenable to an analytic 
treatment that simplifies the numerical computation of 
the forward and backward trajectories, while providing 
sufficient coupling to break the exact intcgrability of the 
harmonic system. 

In the special case k — 0, the system is globally har- 
monic. In this case, the constructions outlined below 
yield a moving dividing surface that is strictly free of re- 
crossings. If fc ^ 0, deviations from the harmonic dynam- 
ics will arise outside the TS region that may lead to error 
in the identification of reactive trajectories. Nonetheless, 
the wealth of microscopic detail that the moving divid- 
ing surface reveals can most easily be illustrated using 
the harmonic limit. This is shown in Sees. ITTTl and FV^ 
The real power of the numerical method, however, lies 
in addressing nonlinear systems; the accuracy of the ap- 
proximate identification of nonlinear reactive trajectories 
is discussed in Sec. IV Bl 

With the potential Q , the Langevin equation ^ reads 



Qa^it) = nq^{t) + 0{ql) - Tq^{t) + ^^{t) 



(5) 



where 



n 







(6) 



is the matrix of second derivatives oiU{q). The nonhnear 
terms in ((S)), which stem from the anharmonic contribu- 
tions to the potential (0)) will be ignored in the remain- 
der of this section, where an exact dividing surface for 
the harmonic limit will be constructed. The full nonlin- 
ear equation of motion (^ will be taken up again in the 
numerical calculations of Sec. IV Bl The following pre- 
sentation can easily be generalized to N spatial dimen- 
sions if y is understood to denote an (n — l)-dimensional 
vector and the corresponding squared frequency ojy an 
(n — l)-dimensional symmetric matrix whose eigenvalues 
need not be degenerate. 



B. The Transition State Trajectory 



As was shown in |29l |30j , Eq. (|5jl can be rewritten in 
phase space, z= {q,v), with v = q, as 



Zait) =AZait) + 







with the 2n-dimensional constant matrix 

/ 

n -r 



(7) 



(8) 



where I is the n x n identity matrix. The matrix A 
is readily diagonalized to yield the eigenvalues ej and 
the corresponding eigenvectors Vj. Equation {Tj) then 
decomposes into a set of 2n independent scalar equations 
of motion 



^ajV') — ^jZaj{t) + SQji^j 



(9) 



where Zaj are the components of ^ in the basis Vj of eigen- 
vectors of A and £,aj are the corresponding components 

of (0,6.(0). 

A particular solution of Eq. is given by 



4w = <! 



ll^e-^^-^^,it + T)dT 

if j such that Re ej < 0, 

if j such that Ree^ > 0. 



(10) 

Whereas a typical trajectory will eventually descend into 
either the reactant or the product wells, the trajectory 
given by Eq. H1(J|) has the important property 29, '^ 
that it remains in the vicinity of the saddle point for all 
time. In this respect it resembles the equilibrium position 
on the saddle that represents the unique trajectory in the 
absence of noise that never descends from the saddle. We 
named this distinguished trajectory the Transition State 



Trajectory in [29|, |30j because it plays as central a role 
in the TST in a noisy environment as the equilibrium 
point docs in conventional TST. Although the integral 
representation (|10|l defines the TS trajectory, it does not 
provide the most efficient way of calculating it. In fact, 
by means of an algorithm that we introduced in [S^l an 
instance of the TS trajectory can be sampled almost as 
efficiently as an instance of the fluctuating force itself. 



C. The relative dynamics 

Once the TS Trajectory 4(i) = (^^W-^KO) is given, 
any other trajectory under the influence of the same noise 
can be described in relative coordinates 



Azit)=(^^fi^=Ut)-ziit) 



(11) 



where the TS Trajectory serves as the origin of a moving 
coordinate system. The relative coordinate vectors can 
be written without a subscript a because they satisfy the 
noiseless equation of motion 



Aq{t) = n A(f(t) - r Ag(t) , 



or, in phase space. 



Az{t) = A Azit) 



(12) 



(13) 



and are, therefore, independent of the noise. Using the 
eigenvectors of A, one can construct invariant manifolds 
and a no-recrossing surface of the noiseless relative dy- 
namics. According to Eq. (|ll|l . they can then be re- 
garded as being attached to the TS Trajectory and being 
carried around by it. In this way one obtains moving in- 
variant manifolds and a moving no-recrossing surface in 
the phase space of the full, noisy dynamics J29ll3(]l |. 

In the two-dimensional case of the potential in Eq. Q) 
under isotropic friction as specified in Eq. ||2Jl, the eigen- 
values of A can be found explicitly: 



eti = -^(7-V7^^4^) , 
et2 = -^ (7+ Y7^-4tJy) • 



(14) 



The corresponding eigenvectors read 



V„ = 



Fti- 



1\ 




/I 







0/ 


Vo 


o\ 


/o 


1 

' 


l^t2 = 
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eJ 




\et2 



(15) 
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FIG. 1: Phase portrait illustrating the relative dynamics 
in the reactive degree of freedom. Dashed lines indicate the 
stable and unstable manifolds of the equilibrium point, and 
the dotted arrows display the corresponding eigenvectors that 
span the diagonal coordinate system. Solid curves illustrate 
representative trajectories. White dots indicate two possibili- 
ties for the position and velocity of the TS trajectory at time 
i = 0, the vertical lines the corresponding barrier ensembles. 
The probability density is given by the line widths. 



These simple analytic results are obtained because iso- 
tropic friction leads to a decoupling of the reactive and 
the transverse degrees of freedom. The eigenvectors Vu 
and Vs span the reactive x-Vx subspace, whereas T4i and 
Vt2 span the transverse subspace y-Vy. 

The knowledge of the eigenvectors allows one to ex- 
phcitly specify the coordinate transformation between 
position- velocity coordinates Ax,Ay,Avx,^Vy and the 
diagonal coordinates Axu, Axg, Axti, Axt2 that charac- 
terize a phase space point via Ax = ^^^AxiVi. In the 
reactive subspace, these transformations read 



Ax — Axu + Accs , 
with the inverse 

AVx - Es Ax 

Ax,, = 



AVx = Eu AXu + Cs AXs (16) 



Ax= 



-AVx + Eu Ax 



(17) 



For all values of 7 and tOx, the eigenvalue eu is pos- 
itive, whereas eg is negative. They correspond to one- 
dimensional stable and unstable subspaces within the re- 
active degree of freedom which, together with represen- 
tative trajectories, are illustrated in Fig. Q] The coor- 
dinate Axu determines the fate of a trajectory in the 
remote future: Trajectories with Aa:u > descend into 
the product well, those with Ax^ < into the reactant 



well. Similarly, a stable coordinate Axg > indicates a 
trajectory that comes out of the product well in the dis- 
tant past, whereas a trajectory with Aoig < comes out 
of the reactant well. A forward-reactive trajectory that 
changes from reactants to products is thus characterized 
by Axg < and Axu > 0, whereas a backward-reactive 
trajectory has Aoig > and Axu < 0. Each reactive 
trajectory crosses the line Ax = once and only once. 
This line, or in several degrees of freedom the hyper- 
surface defined by this condition, can therefore serve as 
a recrossing-free dividing surface between reactants and 
products. Furthermore, the invariant stable and unstable 
subspaces themselves act as separatrices between reactive 
and nonreactive trajectories. Once the initial condition of 
a trajectory is known relative to these separatrices, it can 
unambiguously be classified as reactive or non-reactive. 



III. THE BARRIER ENSEMBLE 

The rate calculation of infrequent events — such as that 
in chemical reactions — can be greatly simplified by sam- 
pling trajectories in the transition state reg ion rather 
than in the reactant region [SJ, y^ |3^ |33, [Sg- The 
transition path sampling technique, for example, focuses 
exclusively on reactive trajectories and therefore miti- 
gates the difficulty of studying high-dimensional systems 
m 113, m m Ha]. Nonetheless, the rate of infrequent 
events has long been known to be described by a flux-flux 
correlation function relative to a fixed dividing surface 
[M 13 li^. (But see Ref. E for a recent enhanced- 
sampling strategy to smooth the potential and thereby 
speed up the calculations.) One difficulty in comput- 
ing the correlation function, however, is the need for the 
simulated trajectories to be evolved for very long times 
simply to determine which trajectories are reactive. We 
will show below that the use of the time-dependent TST 
dividing surface may allow one to resolve that question 
in a more efficient way by identifying the nature of the 
trajectory — viz. reactive or not — at significantly earlier 
evolution times. 

In order to sample the reactive trajectories efficiently, 
it is useful to use initial conditions in which all the parti- 
cles are placed on the fixed TST dividing surface {x = 0) 
at i = 0. That choice guarantees that the trajectories will 
cross the surface at least once, but it does not prevent 
them from recrossing it. Consistent with the Boltzmann 
weighting in the flux- flux-correlation function |3a . l44Ll45l |. 
the initial conditions are distributed along the stable 
transverse direction y and the velocities according to the 
probability density function. 



f{x,y,Vx,Vy 



(2n kBT)^^^ LUy X 
'4) 



exp {-{u^y + vl+ vl)l2 k^T] Six) . (18) 



This choice defines the barrier ensemble. The integration 
of the Boltzmann-weighted flux of these states gives the 
TST estimate of the numerator of the rate expression. 



If all these states were reactive and never recrossed (re- 
turned to) the fixed TST dividing surface, this estimate 
would be exact. The questions to be resolved below con- 
cern the deviation of the true dynamics from this TST 
estimate. These question will be investigated for both 
harmonic and anharmonic barriers. In all cases, the ini- 
tial conditions will be sampled from the same barrier en- 
semble (|TH|) . 

A stochastic trajectory is determined not only by its 
initial condition, but also by the specific instance of the 
fluctuating force that is acting upon it. In a full-fledged 
rate calculation, for example, an average has to be taken 
over both the initial conditions and the noise. The fo- 
cus of this paper, however, is the information that can 
be obtained about the microscopic reaction dynamics us- 
ing the moving TS surface. For simplicitly, a particular 
instance of the noise has therefore been used to illus- 
trate most of the results. Nevertheless, the calculations 
presented here were repeated for several such noise se- 
quences always leading to the same qualitative conclu- 
sions and thereby confirming that the results shown here 
are indeed typical. (These are not shown here for the 
sake of brevity.) While averages over the noise will tend 
to wipe out much of this microscopic detail, it is instruc- 
tive to confirm the convergence of the identification of 
trajectories in calculating averages. In what follows, the 
average of the forward and backward reaction probabil- 
ity will be used to illustrate the convergence and degree 
of accuracy achievable using the moving TS surface to 



identify the reactive trajectories. 



IV. ANALYTIC RESULTS 

Although anharmonicities have to be addessed in a 
typical chemical system, it is helpful to begin with the 
harmonic limiting case because it is amenable to an an- 
alytic treatment. On the one hand, the harmonic limit 
illustrates the level of microscopic detail in which the 
moving TST method allows one to describe the reaction 
mechanism. On the other hand, the analytic results de- 
rived here provide a benchmark against which the per- 
formance of the numerical calculations of Sec. El can be 
assessed. 



A. Reaction probabilities 

The fate of individual trajectories in the barrier ensem- 
ble H18|) can easily be determined if their initial conditions 
are transformed into relative coordinates. The projection 
onto the reactive degree of freedom is illustrated in Fig.^ 
Since in space-fixed coordinates the barrier ensemble is 
centered around q — v — 0, the distribution function in 
relative coordinates is peaked at the stochastic position 
H^q — — (7|(0), Aw = — 1/|(0). It reads exphcitly 



J 



/„i(Aa;, Ay, Aw„ At-^,) = (2^ k^Tf'^ u. 



exp {-{^li^y + vi? + (A«. + vis + i^vy + 4 J')/2 ksT} S{Ax + 4) . (19) 



The forward-reactive part of the ensemble is formed 
by those trajectories whose initial velocity Av^ is so large 
that the trajectory lies above both the stable and the un- 
stable manifold of the equilibrium point. The knowledge 
of the eigenvectors l(T^ allows one to locate these separa- 
trices quantitatively. Reactive trajectories are thus found 
to be characterized by the condition 



Av^ > Avx, 



X^Cs 



a;ji >0 
2;^ <0 



(20) 



Therefore, the probability for a member of the barrier 
ensemble to be forward-reactive is given by 



P{ = / dAx / dAvx X 



/ dAy / dAvy f^ci{Ax,Ay,Avx,Avy 



/ dAv^e-K^{-{Av^+vlS/2kBT} 

1 „ f Av 



1 crfc ( =:^^:]^}^^^:^ 

2 V V2k^ 



(21) 



which has been written in terms of the complementary 
error function l47l 



2 1'°° 

erfc(a;) = —= / exp(-i^) dt 

Vt^ Jx 



(22) 



In a similar manner, backward-reactive trajectories sat- 
isfy 



(23) 




and their probability in the ensemble is 



Ph = - erfc 



At 



'x,raa.x. > ^xa 



(24) 



B. Reaction times 

In contradistinction to a space-fixed dividing surface, 
the moving TS surface is crossed once and only once by 
each reactive trajectory. This allows us to define a unique 
reaction time At^ for each reactive trajectory: It is the 
time when the trajectory crosses the dividing surface, rel- 
ative to the initial time when the coordinates are spec- 
ified by the barrier ensemble. If the initial conditions 
Aa;u(0) and Aa;s(0) in the reactive degree of freedom are 
prescribed, the reaction time can be calculated explicitly. 
The dynamics of the reactive degree of freedom is given 

by 



Ax^{t) = Axii(0)e'"* , Axsit) = A2;s(0)e' 



(25) 



The dividing surface is characterized by the condition 
Aa; — 0, which can be rewritten in relative coordinates 
as Axu — — Axg. The reaction time At^ at which this 
condition is satisfied is easily found to be 



Ati 



1 



1 



In 



In 



-Ax,(0) 
A2;u(0) 

Ai.,(0)-£„Ax(0) 
Ai;:r(0)-esAa;(0) 



(26) 



It is defined for all initial conditions that are either 
forward- or backward-reactive. For a forward-reactive 
trajectory, Az;a;(0) > 0. Because e^ > and eg < 0, it can 
easily be seen from Eq. ^^ that Ai* > if Aa;(0) < 0, 
as it should be for trajectories that start on the reac- 
tant side of the dividing surface and are still to cross 
it. Similarly, a trajectory with Ax(0) > is already on 
the product side, and its reaction time is negative. A 
backward-reactive trajectory, on the other hand, has an 
initial velocity Awj:(0) < 0. In this case, At^ < if 
Aa;(0) < and Ai* > if Aa;(0) > 0. 

If the initial position Aa;(0) is fixed, the reaction 
time (|26|l tends to zero as Awj;(0) ~> oo: Trajectories 
with large initial velocities cross the barrier fast. On 
the other hand, as the separatrices that bound the re- 
active region are approached, i.e. Awa;(0) -^ euAa::(0) if 
Aa;(0) > or Au^(O) -^ esAa;(0) if Aa;(0) < 0, trajecto- 
ries keep barely enough energy to cross the barrier, and 
their reaction times tend to -f oo or — oo, respectively. 

Once the reaction time is given as a function of initial 
conditions, the distribution for the forward- or backward- 
reactive part of the barrier ensemble (|18|l is readily ob- 
tained. In the former case, its probability distribution 



function is given by 

p{At) = — dAx I dAv^ X 

/ dAy / dAvy f^c\{,Ax,Ay,Avx,Avy) x 

5{At-AtHAx,Av.,,Ay,Avy)). (27) 

The normalization factor 1/Pf accounts for the fact 
that only the forward-reactive part of the ensemble con- 
tributes to the distribution. 

The distribution function (|27|l can in its most conve- 
nient form be written in terms of the dimensionless scaled 
time At^ = (eu — es)At-'-. It then reads 



1 



p{At^) =^^p(ArV(eu - Cs)) 



^At' 



r 



exp 



1 -e 



-At* 



where 



w ■ 



gi(0)(eu-es) 



4(0)-eugi(0) 
\/2fc^ 



(28) 



(29) 



(30) 



The reaction probability Pi can be written in terms of r 
and w as 



Pf = 



er 



■fc(r + w) 



^ erfc( 



w) 



r < , 



(31) 



The valid range of Ar* is < At* < oo if q|(0) > 
and -oo < At* < if g*(0) < 0. The distribution 
function H28|l is normalized so that its integral over that 
range is one. Remarkably, the distribution depends only 
on the two parameters r and w, even though the system 
dynamics and the distribution of initial conditions are 
determined by the five parameters Wb, 7, T, q\{Q) and 

In a similar manner, the distribution of reaction times 
can be computed for the backward-reactive part of the 
ensemble. The result is again given by Eq. I|28|) . except 
that the vahd range is now — oo < At* < if q|(0) > 
and < At* < oo if q* (0) < 0. To obtain the proper 
normalization, the reaction probability Pf in the prefac- 
tor of Eq. H28|) must be replaced by the backward-reaction 
probability Ph^ which in terms of the scaled parameters 
reads 



i erfc(— tu) 
5 erfc(- 



w) 



r >0, 
r < 0. 



(32) 
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FIG. 2: The distribution 12811 of reaction times (a) for w = \ 
and r = 0.2 (solid), r — 0.5 (dashed) and r — 1 (dotted), 
(b) for TO = — 1 and r — —0.2 (soUd), r = —0.5 (dashed) and 
r = -1 (dotted). 



As can be seen from Fig. |21 the reaction-time distri- 
bution (|28|l is highly asymmetric around its peak. The 
probabihty distribution function is flat at r = 0, where 
all of its derivatives are zero. For large \t\, it decays 
exponentially like 



p(At*) « < 



-(r+™)^ 



-At* 



V^Pl 



y^Pt 



if 9^(0) >0, Ar^^^+oo, 
iiqi{0)<0, At^^-oo 



(33) 

Because the distribution is so highly asymmetric, the av- 
erage reaction time will be significantly larger than the 
most probable reaction time that is indicated by the max- 
imum of the distribution function. 



NUMERICAL RESULTS 



In what follows, the initial conditions, at t = 0, are cho- 
sen from the distribution in Eq. (|18|1 . All trajectories are 
evolved forward and backward in time to t — ±Tint/'2, 
using the stochastic integration algorithm introduced by 
Ermak and Buckholz ^^ ^J . For the backward propa- 
gation, the integration scheme was modified as described 
in the appendix. In a conventional calculation of the 
exact rate expresssion, reactive trajectories are identi- 
fied according to the positions they attain at the start 
and end of the integration interval: Trajectories that at 
t = —T[^t/2 and t = -|-Ti„t/2 are located on opposite sides 
of the space-fixed dividing surface x = are classified 
as forward- or backward-reactive; others are classified as 
nonreactive. This criterion, however, is only reliable if 
the total integration time Tint is sufficiently large. At 
short times, recrossings of the dividing surface introduce 
unavoidable errors. 

An alternative criterion for the identification of reac- 
tive trajectories is obtained if the space-fixed dividing 
surface is replaced by the moving TS surface described 
above. In the most naive implementation, trajectories 
can be classified as reactive if they are on opposite sides 
of the moving TS surface at t — ±Tinf If the moving- 
TS-surface algorithm is used instead. Tint can be reduced 
by as much as a factor of 2 while still obtaining nearly 
accurate results. In addition, given that the moving TS 
surface is exactly free of recrossings in the harmonic limit 
and approximately so in an anharmonic potential, the in- 
tegration can be stopped as soon as a trajectory crosses 
the moving surface: There is no need to follow the trajec- 
tory further and check for recrossings. Therefore, when 
the moving TS surface is used, the actual integration 
time will on average be much smaller than the nominal 
integration time Tint. 

The reliability of this identification is illustrated be- 
low using the two-dimensional saddle point potential of 
Eq. ^with and without anharmonicity, k. In all of the 
numerical calculations, the units are chosen for simplic- 
ity such that /cbT = 1. The friction is isotropic, with 
7 = 0.2 in these units, and selected so as to be near the 
turnover between the energy- and space-diffusion limited 
regimes. Although most of the calculations assume the 
same fixed noise sequence, averages of the forward and 
backward reaction probabilities over the noise are also 
shown below. In the former, the number of trajecto- 
ries is fixed at A^t — 15 000, which is large enough to 
make statistical errors negligible. In the single-noise cal- 
culations on the two-dimensional harmonic barrier, the 
transverse frequency tOy = 1.5, and the barrier frequency 
is set to ujx — l-O- The latter is reduced to w^: = 0.75 for 
the noise-averaging and in the nonlinear cases in order to 
accentuate the nonlinear coupling. 



As soon as the anharmonicities of the potential in a 
realistic chemical system have to be taken into account, 
the equations of motion can no longer be solved analyti- 
cally, and recourse must be taken to numerical methods. 



A. Harmonic Systems 

A typical reactive trajectory and the TS trajectory in 
the harmonic limit (fc = 0), are shown in Fig. |31 Clearly, 




FIG. 3: The evolution of a member of the ensemble and the 
Transition State Trajectory depicted as the gray and black 
hnes, respectively. The underlying potential is included, and 
the fixed transition state a; = is highlighted by the heavy 
black line. The time-independent projection is shown on the 
base of the figure. The sample trajectory is backward-reactive 
since it is a reactant in the future and product in the past. 
As can be seen, the Transition State Trajectory remains in 
the vicinity of the barrier for all times. 
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FIG. 4: Reaction probabilities calculated using the fixed 
dividing surface displayed as a function of total integration 
time. The fractions of forward-reactive, backward-reactive, 
and nonreactive trajectories are shown as the solid, dashed, 
and dotted lines, respectively. The dash-dotted line repre- 
sents the fraction of trajectories that cross the surface only 
once at Tint = 0. In these simulations, Nt — 15000 trajec- 
tories were integrated, the friction constant 7 — 0.2, and the 
barrier frequency is ui^ = i- 



the space-fixed dividing surface x = 0, in contrast to the 
moving TS surface, is crossed many times. The respective 
percentages of trajectories classified either as reactive and 
nonreactive using the fixed dividing surface are displayed 
as a function of integration time in Fig. 0] Because all 
trajectories start on the dividing surface, at very short 
times, every trajectory is classified as either forvifard- or 
backward-reactive. Subsequent recrossings of the transi- 
tion state result in transient fluctuations of the reaction 
probabilities that slowly approach the true, long-time val- 
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FIG. 5: The number of recrossings (dot-dashed line) of the 
fixed transition state normalized by the total number of tra- 
jectories. The components of the total that resulted in a 
forward-reactive or backward-reactive trajectory are shown 
as the solid and dashed lines, respectively. Recrossings that 
resulted in a nonreactive trajectory are displayed as the dot- 
ted line. The simulation parameters are the same as in Fig.|l] 



ues. Figure 0] also shows the percentage of trajectories 
that are nonreactive as well as those that cross the fixed 
dividing surface only once at Tint = 0. The latter com- 
prise the majority of the reaction events, whereas the 
percentage contributed by reactive trajectories is com- 
paratively small. Nevertheless, the fluctuations in the 
computed reaction probabilities that are caused by re- 
crossings are considerable. 

Because recrossings are crucial to the performance of 
the algorithm, it is instructive to analyze them in more 
detail. Figure |S1 shows the average number of recross- 
ings per trajectory as a function of the total integration 
time. The trivial crossing of the dividing surface that 
all trajectories undergo at t = is not included. The 
number increases monotonically as the trajectories cross 
and recross the transition state. Eventually, it reaches a 
plateau as they leave the barrier region and are lost into 
either the product or reactant states. In addition. Fig. [S] 
decomposes the total number of recrossings into those 
recrossings that occur on trajectories that are found to 
be forward-reactive, backward-reactive or nonreactive at 
the given integration time. Because the classification of a 
particular trajectory can change with increasing integra- 
tion time, these contributions are not nronotonic. Most 
prominently, as the number of nonreactive trajectories 
decreases almost to zero at Tint ~ 2 (see Fig. 0J, the 
contribution of nonreactive trajectories shows the same 
behavior. For large integration times, the largest con- 
tribution to recrossings stems from nonreactive trajecto- 
ries, which are bound to recross the dividing surface at 
least once. In fact, a comparison of Fig. ^ and Fig. [51 
reveals that nonreactive trajectories on average recross 
more than three times before they finally leave the bar- 
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FIG. 6: Reaction probabilities as a function of integration 
time calculated using the moving (solid curves) or the fixed 
(dashed curves) dividing surface. The upper set of curves rep- 
resents forward-reactive probabilities, with the lower set de- 
picting the corresponding back reactions. The dotted curves 
indicate the asymptotic values Pf = 0.3332 and Pb ~ 0.4993 
calculated from Eq. (I31II or 13211 . respectively. The simulation 
parameters are the same as in Fig. 0] 



rier region. Most of the reactive trajectories, by contrast, 
do not recross, and their contribution to the recrossing 
statistics is much smaller. Asymptotically, both forward 
and backward reactive trajectories recross on average ap- 
proximately 0.25 times. 

The dynamics is greatly simplified if the moving TS 
surface is used instead of the fixed one. Reaction prob- 
abilities computed using either surface are compared in 
Fig. Because the trajectories start at a distance from 
the moving TS surface, the corresponding rates are zero 
for short integration times. They then steadily increase 
toward the true long-time probabilities. Since the divid- 
ing surface cannot be recrossed, the asymptotic values 
are approached monotonically. The erratic fluctuations 
of the computed reaction probability that the fixed sur- 
face produces are absent if the moving TS surface is used, 
so that a strict lower bound for the reaction probabil- 
ity is obtained even for very short integration times. In 
quantitative terms, the moving TS surface identifies a 
trajectory as reactive if its reaction time At^ lies within 
the integration interval, so that the finite-time reaction 
probability for a forward reaction is given by 



Pf(Ti„t)=PfxProb<^|Att|< 



Ti, 



(34) 



and a similar expression for the backward-reaction prob- 
ability. The reaction probabilities computed from the 
moving TS surface are therefore determined by the dis- 
tribution (|28|l of reaction times. The convergence toward 
the long-time probability is described by the long-time 
tail H33|l of the reaction-time distribution and is expo- 
nentially fast. Indeed, Fig. El shows that reaction proba- 
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FIG. 7: The reaction probabilities averaged over different in- 
stances of the noise on the harmonic potential for four differ- 
ent values of the number of trajectories (TVt) in the ensemble. 
The solid lines depict the results predicted by Eqs. 1)3 l|l or 
13211 . The dashed and dotted lines are the results obtained 
using the respective fixed or moving dividing surfaces. In the 
harmonic case, these two surfaces provide the same results 
and are indistinguishable. For the case of A't — 1000 the light 
dotted lines display the 95% confidence interval with respect 
to the number of noise sequences sampled (N^). The simula- 
tion parameters are the same as those defined in Fig. |5] 



bilities computed using the moving TS surface converge 
much faster than those obtained from the fixed surface. 
Moreover, in cases such as the current problem, in which 
the separatrices between reactive and non-reactive tra- 
jectories are known exactly, the reaction probabilities Pf 
and Pb can be computed a priori, without having to per- 
form any numerical simulations. The values obtained 
from Eqs. H31() and (|32() are also indicated in Fig. 
They agree precisely with the asymptotic probabilities 
obtained from the simulation. Thus, the moving TS sur- 
face can provide accelerated convergence in the rate for 
finite-time computations for linear problems. 

The analytic reaction probabilities, Eqs. H31|l and H32|l . 
for the harmonic barrier represent the limiting values 
that are obtained for one instance of the noise using a 
large number A't of trajectories. To obtain a macroscop- 
ically observable reaction probablity, one has to average 
these results over a large number N^ of realizations of the 
noise. That average cannot be obtained analytically, but 
it can be easily calculated by a numerical quadrature. It 
provides a useful benchmark for the convergence of the 
computational schemes with respect to A^t and A"^. Fig- 
uredillustrates the forward and backward reaction prob- 
abilities, averaged over N^ realizations of the noise, as a 
function of Aj and for different values of A't- The solid 
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FIG. 8; The distribution of reaction probabilities in the har- 
monic hmit calculated from Eq. ^^ or ||35J from N^^ = 20000 
different instances of the noise. The x-symbols display the re- 
sults for forward-reactive probabilities and the o-symbols are 
for backward-reactive. The simulation parameters are the 
same as those defined in Fig. |U] 



and dashed curves are obtained if reactive trajectories are 
identified through the criteria provided by the fixed and 
the moving TS surfaces, respectively. As expected for a 
symmetric barrier, forward and backward reaction prob- 
abilities converge toward the same limit. Moreover, the 
distributions of forward and backward reaction probabil- 
ities agree, as shown in Fig. |S1 For large Nt, the results 
in Fig. [3 agree with the analytic value displayed as the 
dot-dashed curve in the figure's bottom panel. Therein, 
dotted curves are used to indicate the 95% confidence in- 
terval to further illustrate that the simulations are con- 
verging toward the correct limit as expected. 

The simulation results in Fig. [T] that employ the con- 
ventional criterion for identifying trajectories have been 
computed using the large integration time Tint — 21.5, 
to illustrate the exact results within the error bars of the 
number average. However, it should be clear from Fig.^ 
that the moving-TS-surface criterion often identifies re- 
active trajectories in less than half this time, and once 
so identified a trajectory does not need to be integrated 
further. Given that the calculation of the moving sur- 
face itself — which amounts to the calculation of the TS 
trajectory — takes roughly as much computational effort 
as the integration of an ensemble trajectory, computa- 
tional savings can thus be obtained from the use of the 
moving surface whenever the number Nt of trajectories 
per noise sequence is larger than 2. 



B. Nonlinear Systems 

The true test for the usefulness of the moving transi- 
tion state lies in its ability to identify reactive trajectories 
beyond the linear regime. If nonlinearities are present. 
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FIG. 9: Reaction probabilities as a function of integration 
time calculated using the moving (solid line) or fixed divid- 
ing surface (dashed line) for various values of the coupling 
constant. The time step has been reduced to 8 x 10~® for 
convergence and the barrier frequency changed to ujx ~ 0.75 
to accentuate the nonlinearity. For the case k — 0, the results 
for the reaction probabilities as calculated from Eq. (1311 or 
13211 are included as the dotted lines. 



the relative coordinate Hll|) does not achieve a complete 
separation of the relative motion from the motion of the 
TS Trajectory. Therefore, the moving dividing surface 
will not strictly be free of recrossings. However, if the 
nonlinearities are weak, recrossings can be expected to 
be rare. In these cases, the moving dividing surface will 
be recrossing-free to a useful approximation. Indeed, our 
results indicate that its advantages over a fixed dividing 
surface persist well beyond the harmonic limit. 

We investigate the performance of the moving divid- 
ing surface in the example of the potential (0J, with the 
coupling constant k now taking non-zero values. The re- 
action probabilities for several different values of k are 
displayed in Fig. |5| To accentuate the anharmonicity, 
the barrier frequency was reduced to uJx = 0.75 to al- 
low trajectories to spend more time in the barrier region 
before escaping. For the transverse frequency, the value 
ujy = 1.5 was retained. Evidently, for sufficiently long 
integration times the moving transition state provides 
essentially the same result as the fixed dividing surface 
for all values of the coupling constant up to A: = 0.1. 
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However, the reaction probabilities converge toward the 
long time limit monotonically and much faster than those 
computed with the fixed dividing surface. Therefore, the 
computational advantages that the moving surface offers 
in the harmonic limit persist even in the presence of quite 
substantial nonlinearities. Eventually, of course, the use 
of a moving dividing surface based upon the harmonic ap- 
proximation ceases to be meaningful, as can be seen for 
k = 0.5 and fc = 1. For the specific instance of noise used 
in these calculations, the results obtained from the mov- 
ing and fixed dividing surfaces remain in agreement for 
the backward-reactive trajectories, whereas a substantial 
difference arises for the forward-reactive trajectories. As 
is to be expected of any TST scheme, in these cases the 
moving dividing surface overestimates the reaction prob- 
ability because any trajectory that crosses the surface 
is assumed to be reactive, whereas the possibility of re- 
crossings is neglected. Although not shown, a different 
instance of the noise does not change the trends observed 
in Fig. El 

As in the harmonic limit, the computational advan- 
tages of the moving dividing surface in systems with 
moderate nonlinearities stem from the fact that it is ap- 
proximately free of recrossings. This is illustrated in 
Fig. 1101 The average number of recrossings per trajectory 
of the fixed transition state exhibits similar behavior for 
small to moderate values of the coupling constant. It ap- 
proaches approximately one recrossing per trajectory in 
the long-time limit. In these cases, the number of recross- 
ings of the moving dividing surface is so much smaller 
than the corresponding number for the fixed surface that 
it is not visible in the figure. At larger coupling, the num- 
ber of recrossings of the fixed dividing surface does not 
converge to a finite long-time limit, but instead increases 
linearly with the integration time. The onset of a similar 
behavior occurs at approximately the same value of Tint 
for the moving dividing surface as well. 

This increase in the number of recrossings for both the 
fixed and the moving surfaces is caused by a small per- 
centage of trajectories in the ensemble that never leave 
the TS region for negative times, but rather get trapped 
in an oscillation in the stable transverse degree of free- 
dom y. If the value of y is sufficiently large, the reactive 
degree of freedom x in the potential l@J ceases to be un- 
stable but instead behaves as a harmonic oscillator with 
a (possibly large) effective frequency a)^ = fc y^ — w^ . As 
a result of these fast oscillations in the reactive degree of 
freedom, the dividing surfaces are crossed many times. 
This mechanism has been confirmed by a detailed tra- 
jectory analysis, which for brevity we do not show. It is 
a rather peculiar feature of our model potential due to 
the presence of only one higher order coupling term in 
the potential Q). We would not expect such aberrant 
behavior to arise in a typical system. 

It is clear from Fig. |51 that for moderately strong an- 
harmonicities the moving transition state correctly iden- 
tifies the overall number of reactive trajectories. How- 
ever, that number is a macroscopic observable, and it is 
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FIG. 10: The average number of recrossings of the mov- 
ing transition state (solid lines) and the fixed transition state 
(dashed line) normalized by the total number of trajectories 
for given values of the coupling constant. The simulation pa- 
rameters are the same as those defined in Fig. lO] The values 
for the moving transition state are too small to be seen on 
the same scale in the top four panels. 



not immediately clear whether, on a microscopic level, 
individual reactive trajectories are identified correctly. 
The fraction of trajectories that are identified correctly 
by the moving transition state is displayed in Fig. 1111 
where the "correct" identification for a given trajectory 
has been assumed to be given by the fixed dividing sur- 
face for a sufficiently long integration time. In the cases 
of weak to moderate coupling, the classification obtained 
from the moving dividing surface is correct for all tra- 
jectories, but, as expected, it begins to fail for coupling 
strengths around fc = 0.5. The fact that the identifi- 
cation of the backward trajectories is poorer than that 
for the forward trajectories at large k is not surprising. 
The initial distribution — particles located at the naive 
fixed transition state with forward velocity — disfavors 
backward trajectories which must recross the fixed TS at 
least twice more in order to reach the appropriate bound- 
ary conditions. Nevertheless, Fig. 1111 confirms that the 
favorable behavior of the moving surface that is appar- 
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FIG. 11: The fraction of correctly identified trajectories ac- 
cording to the moving transition state. The correct identifi- 
cation of a trajectory is that defined by the fixed transition 
state at the end of the simulation, Tint = 21. The correctly 
identified backward-reactive or forward-reactive trajectories 
are displayed as the solid and dashed lines, respectively. The 
simulation parameters are the same as those defined in Fig. El 



ent in Fig. O indeed reflects a correct description of the 
underlying microscopic dynamics. 

Figure El displays the noise-averaged reaction proba- 
bilities for the fixed and moving dividing surfaces for a 
couphng of fc = 0.1. The results display the same conver- 
gence behavior in both cases, except that for the moving 
TS surface they are shifted to larger values by roughly 
5%. This small error is due to the small percentage of 
trajectories that recross the moving dividing surface, as 
seen for a particular instance of the noise in Fig.^l Be- 
cause the potential barrier described by Eq.^lis symmet- 
ric even for k ^ 0, the average values of the forward and 
backward reaction probabilities Pf and Pb are equal. The 
simulation results converge rapidly, with respect to both 
A^t and N^, toward their limiting value. These results 
demonstrate that the moving TS surface retains its re- 
liability and its computational advantages for moderate 
values of the anharmonicity upon noise averaging as well 
as for a single instance of the noise. 
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FIG. 12: The reaction probabilities across an anharmonic 
potential (with a coupling of fc = 0.1) are shown as a smooth 
function of the number of different instances of the noise (A'j) 
and a discrete function of the number of trajectories (A^t) 
used to represent the ensemble average. The dashed lines and 
dotted lines result from the use of the fixed or moving dividing 
surfaces in the identification of trajectories, respectively. Note 
that the ordinate ranges over the very narrow interval between 
.37 and .41, and hence the converged exact and approximate 
approaches are nearly equal. The simulation parameters are 
the same as those defined in Fig. j^l 



VI. CONCLUDING REMARKS 

We have recently developed an analytic method for 
constructing a time-dependent stochastic dividing sur- 
face that is strictly free of recrossings i23, |23 ■ In the 
present work, it has been shown that this moving divid- 
ing surface can be used to identify reactive trajectories 
reliably in linear and nonlinear systems: In the harmonic 
limit, the moving dividing surface attached to the TS tra- 
jectory is strictly free of recrossings, while in more general 
(nonseparable) cases it is approximately so. The identifi- 
cation of reactive trajectories using the moving dividing 
surface has been seen in this article to be fairly accurate 
even in the presence of large anharmonic coupling. It can 
be obtained in roughly half the time that is required to 
confirm the nature of a trajectory by numerically evolv- 
ing it to its final state. 

In several of the calculations presented in this article, 
observables have been calculated for a particular instance 
of the noise while averaging over the initial conditions 
of the subsystem. In such restricted averages, the use 
of the moving surface reduces the computational cost of 
the calculation by a factor of two or more. A typical av- 
erage of an observable, however, requires one to include 
multiple instances of the noise. When the average is per- 
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formed using the machinery of the moving TS surface, 
the TS trajectory must be generated for each instance of 
the noise. If only one system trajectory is calculated for 
each noise sequence, the computational effort to calculate 
both the sample trajectory and the moving TS surface is 
about the same as calculating a single (longer) trajec- 
tory. Improved CPU performance can still be obtained if 
one recognizes that the average should be taken by sam- 
pling several trajectories for each noise sequence. Apart 
from the insight into the microscopic reaction dynamics 
that the moving dividing surface offers, it consequently 
also provides computational advantages in the calcula- 
tion of macroscopic observables. Moreover, the compu- 
tation can be readily parallelized because the algorithm 
is embarrasingly parallel when sampling across the tra- 
jectories associated with a given noise sequence. (Indeed, 
although not discussed explicitly in the text, the codes 
have been parallelized across several CPUs with near lin- 
ear scaling.) 

In summary, we envision at least two approaches in 
which the TS-trajectory criterion for reactive trajecto- 
ries will be useful in calculating reaction rates: (i) In 
harmonic (or nearly harmonic) systems, the algorithm 
described here provides a formally exact expression for 
the reaction probability given a noise sequence. This 
term and related averages can be used to substantially 
reduce the required computational time because it lim- 
its the numerical effort to a sampling of the noise, (ii) 
In arbitrary anharmonic systems, the criterion can be 
used to reduce the computational effort to calculate any 
correlation function — such as that in the reaction-rate 
expression — that relies on the correct identification of re- 
active trajectories. The rate expression and other related 
observables that can take advantage of the identification 
of reactive trajectories will be calculated in future work. 
As an illustration, the TS-trajectory criterion was seen 
in this work to converge the forward and backward reac- 
tion probabilities even in a fairly anharmonic case. Thus 
the central result of this work is: the moving dividing 
surface can be used reliably and efficiently to identify 
reactive trajectories. 
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been implemented. The backward integration requires 
special care if one wishes to follow the same stochastic 
trajectory both forward and backward in time. The mod- 
ification of the integration scheme that is necessary to 
this end is described here. 

The forward numerical integrator for Langevin equa- 
tions Ha [43 takes the form 



r{t) + civ{t) + C2a{t) + 5r , (35) 

C3v{t) -f Cia{t) + cr-,a{t + St) + Sv , (36) 



r{t + St) 
v{t + St) 



where a{t) is the acceleration caused by the potential of 
mean force, the Ci are numerical coefficients that depend 
on the time step St and the damping constant 7, and the 
random variables Sr and Sv are sampled from a known 
Gaussian distribution. Time reversal in this algorithm 
can be obtained through a shift in time by ^St so that t 
becomes t — St and t + St becomes t. This replacement 
and a simple reorganization leads to 

r{t-St) = r{t) - civ{t - St) - C2a{t - St) - Sr {il) 

v{t-St) = — \v(t) - C4a(t - St) - C5a(t) - Sv] (38) 
C3 

The backward step (p?7jl cannot be evaluated as it stands 
because the acceleration a(t—St) depends on the position 
r{t — St) that is yet to be determined. To circumvent this 
problem, we substitute Eq. (|38|l into Eq. (|37() to obtain 



r{t - St) = r{t) - — [v{t) - c^a{t) ~ Sv] 

C3 



Sr 



C1C4 
C3 



C2 a{t - St) . (39) 



When the acceleration a{t — St) is expressed in terms of 
the position r{t — St) through the equation of motion, 
Eq. 139|1 becomes an implicit equation for the positions 
r{t — St) at the earlier time. For all but the simplest 
potentials, it cannot be solved explicitly. Specifically, for 
the anharmonic potential Q), 

U{x,y) ^ --uj^x +^^yy +kxy , 
it leads to the coupled equation system 
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APPENDIX: THE NUMERICAL INTEGRATOR 

The simulations of the reaction dynamics presented in 
Sec. El require one to follow a stochastic trajectory nu- 
merically from t = both forward in time to i = T/2 and 
backward in time to t = —T/2. For the forwar d prop aga- 
tion, a standard stochastic integration scheme p8ll49l | has 



x{t~St)^X{t)+ (^1^-C2 I X 

{tulxit - St) ~ 2kx{t - St)y{t ~ St)^) , 

(40) 

y[t-St)^Y{t)-l^-^-c^^ 

{u;ly{t - St) + 2kx{t - St)^y{t - St)) , 

(41) 

where X{t) and Y{t) denote the contributions of the first 
three terms in Eq. H39fl . In the harmonic limit fc = 0, 
the two equations uncouple and can be solved for simple 
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explicit expressions for the position updates. For nonzero 
k, Eqs. l|in|l and iflTl) represent an implicit integration 
scheme. It can be converted into an explicit method by 
rearranging the terms into 



x{t - St) 
y{t - St) 



x{t) 



l-{'-^~C2]{u^l-2ky{t-5tY) 



(42) 



Y{t) 



l+[^^C2){ujl + 2kx{t-5ty 



■(.43) 



The denominators in Eqs. H42|l and (|43|l are updated 
using Eqs. (|40l) and (|41l) . but the unknown corrections 



involving a{t — 5t) are neglected because the coefficient 
C1C4/C3 — C2 is of second order in the time step 5t. This 
leads to 



x{t~ St) 
y{t ~ St) 



X{t) 



1 - ("if - ^2) i^l - ^kY{tr) 
Y{t) 



C1C4 

C3 



C2){u;l + 2kX{ty) 



, (44) 
.(45) 



Finally, we insert these approximations into the right- 
hand sides of Eq. H42|l and (|43|l to obtain an explicit 
integration scheme backwards in time. 
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